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Abstract 

In  this  paper  we  present  a  coordinate-split  (CS)  technique  for  the  numerical 
solution  of  the  equations  of  motion  of  constrained  multibody  dynamic  systems. 

We  show  how  the  coordinate-split  technique  can  be  implemented  within  the 
context  of  commonly  used  solution  methods,  for  increased  efficiency  and  relia¬ 
bility. 

A  particularly  challenging  problem  for  multibody  dynamics  is  the  numeri¬ 
cal  solution  of  highly  oscillatory  nonlinear  mechanical  systems.  Highly  stable 
implicit  integration  methods  with  large  stepsizes  can  be  used  to  damp  the  os¬ 
cillation.  if  it  is  of  small  amplitude.  However,  the  standard  Newton  iteration 
is  known  to  experience  severe  convergence  difficulties  which  force  a  restriction 
of  the  stepsize.  We  introduce  a  modified  coordinate-split  (CM)  iteration  which 
overcomes  these  problems.  Convergence  analysis  explains  the  improved  conver¬ 
gence  for  nonlinear  oscillatory  systems,  and  numerical  experiments  illustrate  the 
effectiveness  of  the  new  method. 
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1  Introduction 


The  equations  of  motion  of  a  constrained  multibody  system  can  be  written  as  [8] 

q  —  v  =  0  (1.1a) 

M(q)v  +  GT\  -  f(v,q,t)  =  0  (1-lb) 

g(q)  =  0  (1-lc) 

where  q  =  [q\,  q 2,  qn]  are  the  generalized  coordinates,  A  =  [Ai,  A2, Am]  are  the 
Lagrange  multipliers,  M(q)  €  RnXn  is  the  mass-inertia  matrix,  /  G  2R”  is  the  force 
applied  to  the  system,  q  =  is  the  velocity  and  q  =  ^  is  the  acceleration  vector. 
The  constraints  g  =  [pi,  pm]  are  m  smooth  functions  of  q,  whose  Jacobian 

G(q)  =  NN  €  5?mXn  ,  m  <  n  (1.2) 

|dPj. 

is  assumed  of  full  row-rank.  We  assume  that  G{q)M{q)GT(q)  is  symmetric  and  posi¬ 
tive  definite  for  every  q  €  lRn  to  obtain  a  consistent  physics  represented  by  (1.1).  The 
degrees  of  freedom  for  the  system  (1.1)  is  n  —  m.  Equation  (1.1)  is  a  well-known 
index-3  DAE  [3,  11]. 

Many  methods  have  been  proposed  for  modeling  multibody  systems.  Direct  nu¬ 
merical  integration  of  the  index-3  DAE  (1.1)  suffers  from  the  well-known  difficulties 
inherent  in  the  solution  of  high-index  DAEs  [11].  One  way  to  lower  the  index  involves 
introducing  derivatives  of  the  constraint  g{q),  along  with  additional  Lagrange  mul¬ 
tipliers  p.  This  yields  the  stabilized  index-2  or  GGL  formulation  of  the  constrained 
equations  of  motion  [6] 


q  —  v  +  GTp  =  0 

(1.3a) 

M(q)v  +  Gt\  —  f(v,  q,  t)  =  0 

(1.3b) 

O 

II 

(1.3c) 

ST 

II 

0 

(1.3d) 

which  has  been  widely  used  in  simulation.  The  Lagrange  multiplier  variables  A  and 
p  fulfill  the  role  of  projecting  the  solution  onto  the  position  (1.3d)  and  the  velocity 
(1.3c)  constraints,  respectively.  Equations  (1.3)  and  related  systems  have  been  solved 
by  a  variety  of  methods.  Here  we  will  consider  solution  by  implicit  numerical  methods 
such  as  BDF  or  RADAU.  A  closely  related  approach  is  based  on  explicitly  projecting 
the  numerical  solution  onto  the  constraints  [16,  18,  20,  21]  and  involves  many  of  the 
same  issues  for  the  implementation  that  are  considered  here. 

Many  of  the  numerical  methods  for  multibody  systems  solve  the  system  (1.3) 
directly.  It  is  also  possible  to  eliminate  the  Lagrange  multipliers  and  reduce  the  size 
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of  the  system  to  the  number  of  degrees  of  freedom.  One  way  to  accomplish  this 
begins  with  the  stabilized  index-2  system  (1.3).  Suppose  that  G(p)  is  full-rank  on  the 
constraint  manifold  M  =  {q  (E  FLn  \  g(q )  =  0}.  Then  one  can  find  an  annihilation 
matrix  P(q )  €  J?(n-m)Xn  such  that  P(q)GT(q)  =  0,  €  M.  Premultiplying  (1.3a) 

and  (1.3b)  by  P(q)  yields  an  index- 1  DAE 


P(q)(q  -  V) 

=  0 

(1.4a) 

P(q)(M(q)v  -  f(v,q,t)) 

=  0 

(1.4b) 

G(q)v 

=  0 

(1.4c) 

9{q) 

=  0. 

(1.4d) 

There  is  a  potential  gain  in  efficiency  for  this  formulation  due  to  the  size-reduction 
of  the  nonlinear  system,  compared  to  (1.3).  An  important  practical  consequence  of 
(1-4)  is  that  (p,  A)  have  been  eliminated  from  the  DAE,  via  multiplication  of  (1.3a, 
1.3b)  by  the  nonlinear  P(q).  Thus,  the  error  test  and  Newton  iteration  convergence 
test  in  a  numerical  implementation  of  (1.4)  no  longer  need  to  include  (p,  A).  These 
higher-index  variables  can  cause  problems  in  the  direct  numerical  solution  of  (1.3). 
One  could  in  principle  also  consider  removing  (p,  A)  from  the  test  in  the  solution 
of  (1.3),  however  it  is  not  usually  possible  to  justify  this  action,  particularly  in  the 
case  of  the  Newton  convergence  test.  Elimination  of  these  variables  from  the  Newton 
convergence  test  in  the  solution  of  (1.3)  can  lead  to  a  code  which  sometimes  produces 
incorrect  solutions.  It  is  the  fact  that  multiplying  by  the  nonlinear  P(q)  eliminates 
(p.  A)  from  the  nonlinear  system,  which  allows  these  variables  to  be  excluded  from 
the  tests  in  the  solution  of  (1.4). 

Direct  numerical  solution  of  (1.4)  presents  some  challenges.  First  we  must  have 
a  means  of  generating  P(q)  which  is  reliable  and  cheap.  Further,  we  note  that  the 
Jacobian  matrix  for  the  Newton  iteration  involves  complicated  terms  which  arise 
from  the  derivatives  of  P(q).  We  need  a  means  of  generating  the  Jacobian  matrix. 
Finally,  practical  issues  such  as  the  error  test  and  Newton  convergence  test  must  be 
considered. 

In  the  first  part  of  this  paper,  we  show  how  the  numerical  solution  of  (1.4)  can 
be  accomplished  reliably  and  efficiently.  In  Section  2,  we  show  how  to  obtain  a 
cheap  representation  for  P(q ),  and  then  how  to  compute  the  Jacobian  matrix  without 
directly  computing  the  complicated  derivatives  of  P.  We  show  that  the  nonlinear 
iteration  converges.  The  effectiveness  of  this  method  for  mechanical  systems  will  be 
demonstrated  in  numerical  experiments  in  Section  5. 

Our  approach  for  obtaining  a  cheap  representation  of  P(q)  is  based  on  a  coordinate¬ 
splitting  of  the  variables.  A  widely-used  method  which  is  related  in  the  sense  of  also 
making  use  of  a  splitting  of  the  coordinates  is  the  generalized  coordinate  partitioning 
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(1.5) 


method  [20],  which  yields  n-m  differential  equations 

M(x,h(x))x  =  f(x,^-x,x,h(x),t) 
ax 

where  q  =  Xx  +  Yy  such  that  X  €  2Rnxp  and  Y  €  Rnxm,  whose  columns  constitute 
the  standard  basis  for  2Rn.  The  matrix  Y  is  selected  so  that  ( G(q)Y)~ 1  exists  in  a 
neighborhood  of  q ,  and  h(x)  is  the  implicit  function  of  y  defined  by  the  constraints. 
However,  this  differs  substantially  from  the  approach  we  outline  here  because  P(q) 
associated  with  this  method  is  not  orthogonal  to  Cffiq).  Hence  the  index-reduction  by 
differentiating  the  constraints  and  projecting  to  the  invariant  space  must  be  carried  out 
explicitly.  In  particular,  this  requires  forming  the  derivative  of  the  velocity  constraints 
(i.e.,  the  acceleration  constraints)  explicitly.  Another  method  for  (1.4)  has  been 
proposed  by  [7,  16,  17,  18],  where  V  is  chosen  to  be  an  orthonormal  basis  of  the  local 
tangent  space  of  the  constraint  manifold.  Choosing  a  smoothly  varying  V  is  required 
and  may  cause  some  practical  difficulties. 

Direct  numerical  solution  of  (1.4)  via  our  coordinate-split  approach  yields  an  effi¬ 
cient  and  reliable  method  for  solving  equations  of  motion  for  most  multibody  mechan¬ 
ical  systems.  However,  there  is  a  class  of  multibody  systems  which  present  additional 
computational  challenges.  These  are  the  problems  with  high-frequency  nonlinear  os¬ 
cillations.  Highly  oscillatory  components  are  often  used  to  model  devices  with  strong 
potential  energy.  Typical  examples  of  such  problems  arise  from  modeling  flexible 
multibody  mechanical,  and  molecular  dynamic  systems.  For  many  problems,  oscil¬ 
lations  of  a  sufficiently  small  amplitude  are  not  important  for  the  model,  but  they 
severely  restrict  the  stepsize  for  numerical  methods.  For  these  types  of  problems, 
stiffly  stable  implicit  numerical  integration  methods  can  be  used  to  damp  out  the  os¬ 
cillation  [15].  However,  the  stepsize  may  still  be  severely  restricted  due  to  difficulties 
in  converging  the  Newton  iteration  for  larger  stepsizes  [15].  We  have  studied  this 
class  of  oscillating  problems  in  [22].  The  solutions  are  composed  of  a  low-amplitude 
high-frequency  oscillation  around  a  smooth  solution  [19].  Along  the  smooth  solution, 
the  eigenstructure  of  the  local  Jacobian  matrix  varies  smoothly.  However,  along  the 
solutions  which  are  nearby  to  the  smooth  solution,  the  local  eigenstructure  oscillates 
with  the  high  frequency,  and  is  very  badly  behaved.  The  standard  Newton  iteration 
inside  a  damping  numerical  method  starts  from  a  predictor  which  is  on  a  nearby  so¬ 
lution,  and  attempts  to  find  the  smooth  solution.  It  evaluates  its  Jacobian  matrix 
on  the  nearby  solution,  which  determines  the  direction  it  takes  toward  the  smooth 
solution.  Unfortunately,  these  Jacobian  matrices  do  not  yield  good  directions  for 
nonlinear  oscillating  problems  as  described  above,  unless  the  predictor  is  already  ex¬ 
tremely  close  to  the  smooth  solution.  Thus,  the  standard  Newton  method  must  be 
coupled  with  a  severe  reduction  in  the  timestep  to  achieve  an  adequate  predictor. 

In  Section  3  we  introduce  a  modification  to  the  Newton  iteration  which  we  call 
the  CM-iteration.  This  iteration  is  easy  to  implement,  effective  for  non-oscillatory 
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problems,  and  particularly  effective  for  nonlinear  highly  oscillatory  problems.  The 
basic  idea  of  the  CM-iteration  is  that  there  are  terms  in  the  Jacobian  which  involve 
derivatives  of  the  projection  onto  the  smooth  solution.  These  terms  are  complicated 
to  compute,  large  and  oscillatory  away  from  the  smooth  solution,  and  zero  on  the 
constraint  manifold.  The  CM-iteration  sets  these  terms  to  zero,  yielding  a  reliable 
direction  towards  the  smooth  solution  for  the  Newton-type  iteration.  In  Section  3, 
we  outline  the  CM-iteration,  give  details  for  its  implementation,  and  prove  its  conver¬ 
gence.  In  Section  4  we  describe  in  more  detail  the  structure  of  nonlinear  oscillatory 
mechanical  systems,  and  derive  estimates  for  the  rates  of  convergence  of  the  CS  and 
CM-iterations  applied  to  these  oscillatory  systems.  The  difference  in  convergence 
rate  explains  why  the  CM-iteration  is  highly  effective  for  oscillatory  systems,  and 
shows  that  its  rate  of  convergence  for  non-oscillatory  systems  is  similar  to  that  of 
the  CS  iteration.  In  Section  5,  numerical  experiments  are  given  which  demonstrate 
the  effectiveness  of  these  methods,  particularly  for  oscillatory  nonlinear  mechanical 
systems. 


2  The  Coordinate-Split  Technique 


In  this  section  we  present  the  coordinate-split  (CS)  technique.  We  show  how  to  define 
the  matrix  P(q )  via  a  coordinate  splitting  and  how  to  compute  this  matrix  cheaply. 
Although  at  first  glance  it  would  appear  that  implementation  of  this  method  would 
be  difficult  due  to  complications  in  computing  the  Jacobian  of  P(q),  we  show  that  the 
special  form  of  the  pseudo-inverse  can  be  used  to  give  a  much  simpler  derivation  of  the 
Jacobian.  We  outline  an  efficient  implementation  for  solving  the  nonlinear  system, 
and  observe  that  the  CS  technique  leads  to  a  natural  and  effective  error  estimator. 
Finally,  we  prove  the  convergence  of  the  CS  iteration. 


2.1  Construction  of  P(q) 

The  construction  of  the  annihilation  matrix  P(q)  involves  the  solution  of  a  class 
of  pseudo-inverses  of  the  constraint  Jacobian  G(q).  An  effective  way  to  obtain  the 
projected  vector  P(q)r  is  to  use  a  splitting  of  the  original  coordinates. 

Definition  2.1  [Coordinate-Splitting  Matrix]  Let  X  and  Y  be  the  matrices  whose 
columns  constitute  the  standard  basis  of  Mnxn  such  that  ||(G(5)Y)-1 1|  is  bounded  in  a 
neighborhood  Uq  of  qo.  The  pxn  coordinate-splitting  matrix  for  (1.1)  is  defined  by 

P(q)  =  XT-  Q(q)TYT  =  XT(I  -  G(q)T(G(q)Y)-TYT )  (2.1) 

where  Q(q)  =  (G(q)Y)-‘G(q)X. 
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Remark  2.1  From  the  construction  of  the  CS  matrix  P(q),  one  can  easily  see  that 
P(q)GT(q )  =  0  for  all  q  €  TRn,  ie.,  P(q)  is  orthogonal  to  range(GT).  Furthermore, 
the  row  vectors  of  P(q)  are  orthonormal,  i.e.,  P(q)TP(q )  =  Ip  where  Ip  is  the  identity 
matrix  in  JRP. 


The  computation  of  P(q)  can  be  carried  out  using  the  LU-factorization  of  the 
constraint  Jacobian  matrix.  Applying  Gaussian  elimination  with  row-pivoting  to  GT 
yields 

£m--.£1GT  =Lm---L,U  (2.2) 

where  Et  is  the  elementary  permutation  and  Li  is  a  Gauss  transformation,  i  € 
(1, 2, m}.  From  the  factorization,  we  have 

[Y,X}  =  E  =  Em---E1.  (2.3) 

Using  the  standard  solution  technique  by  LU-decomposition,  the  projected  vector 
P(q)r  can  computed  in  a  straightforward  and  relatively  cheap  way.  In  addition, 
the  derivative  (2.9)  can  be  computed  using  the  same  factorization  of  GT  and  the 
intermediate  result  s  =  ~(GY)~1YTr  from  the  computation  of  P(q)r. 


Remark  2.2  Alternatively,  one  can  apply  QR-, factorization  to  GT  for  the  computa¬ 
tion  of  P(q)r.  Using  QR- factorization  with  partial  column  pivoting  [10],  we  obtain 

GtE \  •  •  •  Em  =  U  •  •  •  LmU  (2.4) 

where  Ei  is  the  elementary  permutation,  Li  is  the  Householder  matrix,  i  €  {l,2,...,m}, 
and  U  is  upper  triangular.  The  last  (n  —  m )  columns  of  the  orthogonal  matrix 
L  =  nr=,  Ti  constitute  a  basis  for  the  null-space  of  G.  Thus  we  can  write 

[Y,X]  =  £  =  L1.-.Lm.  (2.5) 

Note  that  X  and  Y  are  usually  subsets  of  the  standard  basis  in  lRn . 


2.2  Computation  of  Jacobian  and  projected  vector 

The  derivative  of  (2.1)  can  be  obtained  from  the  formulas  given  in  [9]  (Theorem  4.3, 
pp.  420).  Here,  we  give  a  much  simpler  derivation  of  the  Jacobian  of  P(q),  because 
of  the  special  form  of  the  pseudo-inverse  in  P(q).  For  G(q)  a  smooth  function,  we 
can  in  addition  derive  an  approximation  of  the  projected  vector  function  P{q)r  using 
a  nearby  point. 
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Lemma  2.1  Suppose  A(z)  €  ]RNxN  is  a  nonsingular  matrix-valued  function  of  z  € 
2R‘V ,  whose  components  are  non-constant  smooth  functions  of  z.  Then  the  derivative 
of  A{z)~1r  may  be  obtained  as 

fj  A  ]■  (  rl 

- — - =  —A~1(z)—[A(z)w],  with  w  =  A~lr  (2.6) 

dz  dz 

where  r,w  €  JRN . 


Proof.  Let  r(r)  be  a  smooth  vector- valued  function.  Then  the  derivative  of  A  1(z)r(z) 
with  respect  to  z  leads  to 


f(A-'(z)r(z))  =  ±A->(z)r  +  A-\z)^ 

according  to  the  product  rule.  Choosing  r(z)  =  A(z)w  yields  A~1(z)r(z) 
constant  vector,  and  (2.7)  becomes 


where  w  is  A  lr.  □ 


(2.7) 
=  w  a 


(2.8) 


Lemma  2.2  The  derivative  of  P(q)r  is 

—P(q)r  =  P(p)  ~ ,  with  s  =  -(GY)~TYTr  (2.9) 

where  P(q)  is  defined  by  (2.1)  and  r  6  lRnXl. 


Proof.  Differentiating  P(q)r  with  respect  to  q  yields 

fqP(l) r  =  fqXT(l  -  QT(q)YT)r.  (2.10) 

Since  vanishes,  differentiating  (2.10)  by  the  product  rule  yields 

where  s  =  —  ( GY)~TYTr .  According  to  (2.6)  the  second  term  on  the  right-hand  side 
of  (2.11)  becomes 


{gx)t&^X1l  =  qt*{GY)t* 

dq  dq 

This  may  be  substituted  back  into  (2.10)  to  obtain  (2.9).  □ 
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2.3  Solving  the  nonlinear  system 


Here  we  outline  an  efficient  implementation  of  the  coordinate-split  iteration  for  solving 
the  nonlinear  system  at  each  time  step,  and  show  how  the  coordinate  splitting  leads 
to  a  natural  and  reliable  error  estimator.  In  particular,  the  local  error  estimators 
based  on  differences,  i.e.,  the  increments  of  the  nonlinear  iterations  for  the  discretized 
nonlinear  equations,  can  be  obtained  from  those  on  the  independent  variables  only. 

Applying,  for  example,  a  BDF  formula  to  (1.4)  yields  the  nonlinear  system 

PM(Pk<l«  -On)  =  0  (2.12a) 

PM(MMPhv«  -/(»„,  =  0  (2.12b) 

<?(?„)  o,  =  0  (2.12c) 

9(9.)  =  0  (2.12d) 

where  ph  is  the  discretization  operator,  and  h  the  stepsize  of  the  time  discretiza¬ 
tion.  Given  an  initial  prediction  applying  Newton-type  methods  to  (2.12) 

requires  the  solution  of  linear  equations 


^(^Vn^n  )(^9ni^^’n)  —  ^{Qni^n) 

such  that  A qn  and  Avn  are  the  increments  of  qn  and  vn, 


(2.13) 


?(<?n)TSl  .  dph  Qn] 

dqn  dqn  J 


j ^n)  — 


dqn 

GM 


-PM 

PM— 

G(qn) 

o 


(2.14) 


^(9n?^n)  —  [F* n^l  t  Pn^2i  GnVn,  Qn\ 

where  sx  =  -(GY)~TYTr1,  s2  =  ~(GY)~TYTr2 ,  n  =  phqn~v „,  and  r2  =  M{qn)phv„- 
f(vn,qn,tn).  The  subscript  n  of  a  function  represents  its  numerical  value  at  tn ,  e.g., 

9n  =  p(9n)- 

To  solve  the  Newton  equations  (2.13)  efficiently,  we  apply  a  decoupling  technique 
to  the  equations.  To  begin,  we  rewrite  the  first  two  equations  of  (2.13),  i.e.,  corre¬ 
sponding  to  the  derivatives  of  (2.12a)  and  (2.12b), 

[o{,n)  [a?; ]  +  [!!])  =°  (215) 


where  the  2 n  x  2 n  matrix  J\  is 


{Qm  ^n) 


dPh1«  4.  -J 

dqn  dqn 

i^)PhVn)  _  d[n  i  d(GnS*)  Mdf>hVn  -  ^ 
dqn  dqn  ~  dqn  dvn  dvn 


(2.16) 
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Since  J\  is  generally  invertible  under  the  assumption  of  M(qn)  nonsingular,  one  so¬ 
lution  of  (2.15)  can  be  computed  by 


(XT-QlYT)Aq„  =  0  (2.17a) 

(XT-QlYT)  A0„  =  0  (2.17b) 


where 


Vn) 


A  qn  1 

’  -n  ' 

Avn  _ 

.  r2 . 

(2.18) 


Note  that  the  solution  of  (2.17a)  is  not  the  necessary  but  only  a  sufficient  one  to 
(2.15),  i.e.,  there  are  other  solutions  to  (2.15).  Nevertheless,  we  will  show  in  what 
follows  that  (2.17a)  is  consistent  with  the  iterations  for  the  constraints  along  the  same 
direction. 


Observing  that  A qn  can  be  determined  by  (2.17a)  and  the  last  row  of  (2.13), 
i.e.,  GnAqn  =  -gn,  the  resulting  increments  can  be  used  to  compute  ^( qn)Aqni 
such  that  Avn  is  uniquely  determined  by  (2.17b)  and  the  third  row  of  (2.13),  e.g., 
GnAvn  =  —Gnvn  —  “j^( qn)Aqn .  Thus,  (Aqn,  Avn)  can  be  obtained  by  solving  two 
linear  systems  of  the  form 


Pn  1 

’  a  ' 

.  Gn 

u  = 

b 

for  u  €  JRn,  a  €  and  b  €  lRm.  It  is  important  to  note  that  (2.19)  represents 

Gauss-Newton  iteration  along  the  column  space  of  PT.  Since  the  a- vector  in  the  right- 
hand  side  of  (2.19)  has  been  computed,  we  can  then  compute  u  with  b2  =  PnAqn  for 
Aqn  and  b2  =  Pn Avn  for  Avn. 

Denoting  u  =  Xux  +  Yuy  and  p  =  n  —  m,  the  first  p  equations  of  (2.19)  are  in  the 
form 

Q  Uy  —  o,  (2.20) 

and  the  second  m  equations  yield 

uy  =  (GY)-\b-(GX)ux).  (2.21) 

Substituting  (2.21)  into  (2.20)  and  solving  for  ux  yields 

(IP  +  QtQ)ux  =  (a  +  QT(GY)~1b)  (2.22) 

where  Ip  is  the  p  x  p  identity  matrix.  Using  the  consistent  projection  vectors  PnAqn 
and  PnAvn ,  we  obtain 

X  &qn  =  Axn  =  —  (Ip  +  QnQn)  1(PnX 9n  +  Q^(Gny)-1^n)  (2.23a) 

YTAqn  =  Ayn  =  -(GnY)~l(gn  +  (GnX)Axn)  (2.23b) 
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and 


XTAvn  =  Awn  =  —(Ip+Q^Qn)  ^(PnAvn  +  Q^(GnY)  1rjn)  (2.24a) 
YTAvn  =  Azn  =  -(Gny)-1(r?n  +  (GnX)AW;n)  (2.24b) 

where  gn  =  ^^Aqn  +  Gnvn)  and  vn  =  Xwn+Yzn.  The  numerical  solutions 
(2.23)  and  (2.24)  illustrate  that  the  dependent  variables  y„  and  zn  are  determined 
geometrically  along  the  orthogonal  complement  of  the  column  space  of  Pj.  The 
system  dynamics  is  described  by  the  local  independent  generalized  coordinates  x 
and  its  velocity  w. 

Implementation  of  simplified  Newton  method 

Simplified  New- ton  iterations  are  commonly  used  for  the  iterative  solutions  of  (2.12). 
The  increments  by  Newton’s  method  to  (2.12)  are  (2.23)  and  (2.24).  For  simpli¬ 
fied  Newton  iteration  that  uses  a  fixed  approximation  J(n0)  of  the  Jacobian  J{qn,v „) 
in  (2.13),  the  numerical  solutions  of  (2.23)  and  (2.24)  are  computed  by  some  fixed 
matrices.  For  instance,  (2.23a)  becomes 

Axn  =  -(/,  +  Q<£)TQW)-'(PnAqn  4-  Q(°)r(Gf  y)-1^) 

where  the  superscript  (0)  indicates  the  function  approximated  at  some  (gf0),ui0))- 
Since  (IP+Q^T Q^)  is  symmetric  positive-definite,  the  increments,  under  the  weighted- 
norm  induced  by  diag{Ip  +  Q ^  Im},  become 


Axf  =  PnAfn  +  r(GS5»y)-19« 

A#f  =  -(G^Y)-'(g,  +  (G<„0>X)Axf ) 

A  ws„  =  PnAvs,  +  ^T(G(°)Y)-'4 
Azf  =  -(G<0>y)-1(ih  +  (GL°>Jf)Awf). 


where  Agf  and  A  if  result  from 


(2.25a) 

(2.25b) 

(2.25c) 

(2.25d) 


(2.26) 


Efficient  implementation  of  a  simplified  Newton  iteration  for  (2.12)  is  carried  out 
using  (2.2)  as  follows 

Algorithm  2.1  [CS  iteration] 

Step  1.  Apply  (2.2)  to  GT(q^)  and  Ji{q^\v(^y)  as  in  (2.16). 
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Step  2.  Solve  for  Aq^  and  Av^  in  (2.26). 

Step  3.  Apply  (2.2)  to  GT(qn),  then  compute  PnAq f  and  PnAv„. 

Step  4.  Compute  (G^y)-1^,  where  62  =  0,  then  apply  to  the  result  to  obtain 
(2.25a). 

Step  5.  Use  (2.25a)  to  compute  A t/f ,  where  62  =  Axf . 

Step  6.  Use  (2.25a)  and  (2.25b)  to  compute  r\ then  repeat  Step  4-5  for  (2.25c) 
and  (2.25d). 

For  a  straightforward  implementation,  the  iterative  solutions  use  (2.14)  directly,  which 
requires  the  solution  of  4n  projected  vectors  for  the  first  2(n  —  m)  rows.  For  Algo¬ 
rithm  2.1,  there  is  no  need  to  compute  the  4n  vectors.  Instead,  two  projected  vectors 
and  four  solutions  of  a  m  x  n  linear  system  are  required  for  each  iteration  in  addition 
to  the  solution  of  a  2n  x  2n  linear  system,  e.g.,  (2.14)  and  (2.16)  for  the  straightfor¬ 
ward  and  proposed  implementations,  respectively.  Therefore,  the  above  algorithm  is 
preferred  when  n  is  large.  Note  that  the  computational  cost  of  the  iterative  solution 
by  Algorithm  2.1  is  almost  equal  to  that  of  the  solution  of  a  2(n  +  m)  x  2(n  +  m) 
linear  system.  Using  Lt/-decomposition,  the  solution  of  a  2(n  +  m)  x  2 (n  +  m)  linear 
system  requires  0(4(n  +  m)2)  flops ,  while  the  proposed  iterative  solution  requires 
0(6(nm)2)  +  0(4n2)  flops  with  an  additional  factorization  and  two  solutions  of  the 
matrix  (GnY)~TYT  in  Step  3  in  the  above  algorithm. 

Applying  numerical  integration,  convergence  of  the  iterative  solutions  by  Algo¬ 
rithm  2.1  can  be  achieved  if  the  initial  guess  <3^  from  the  predictor  is  close  enough 
to  the  numerical  solution  qn.  In  addition,  the  local  error  estimator  based  on  the 
predictor-corrector  difference  can  be  modified  for  a  more  effective  approximation. 
From  (2.25b,  2.25d),  the  increments  of  yn  and  zn  are  bounded  by 

IIA^II  <  ||(G<°>y)-,||(||G<0)X||||:r<‘>ll  +  ll<ai) 

IIA4f»ll  <  ii(G<n0>y)-'||(iiG<,”>xiiii»<i>ii  +  u-ai) 

for  all  i.  Ideally,  the  error  estimator  should  be  based  on  the  dynamics  of  such  systems, 
which  is  described  by  x  and  w.  Under  the  assumption  that  the  constraints  are  smooth, 
the  local  error  can  then  be  approximated  using  the  alternative  sequences  of  A and 
A by  setting  <?£)  and  rf)'’  to  zero  in  (2.25b)  and  (2.25d),  respectively.  Since  (gn,  vn) 
converges  in  N  iterations,  the  difference  can  be  estimated  by 

ll(fc,«w)  -  (9'0U0,)II  =  BlIAai  +  l|A»«||)  *  «.E(IIA4?II  +  l|A«4‘>||)  (2.27) 

t=l  *=1 

where  Kn  =  ||(G^y)-1G(l°)A’||.  Therefore,  the  local  errors  of  numerical  integration 
of  (2.12)  can  be  approximated  using  only  the  increments  of  and  tc^,  e.g.,  the 
right-hand  side  of  (2.27). 
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2.4  Convergence  of  the  CS  iteration 


For  simplicity  we  now  consider,  instead  of  the  second-order  constrained  equations  of 
motion  (1.1),  a  first-order  system 

q-f(q,t)  +  GT\  =  0  (2.28a) 

g(q)  =  0  (2.28b) 

since  the  convergence  of  (2.28)  can  be  trivially  extended  to  (1.1).  Applying  stiffly 
stable  numerical  methods,  the  convergence  result  is  well-known,  see  [11]  pp.  494-498. 
Convergence  of  discretization  methods  for  the  index-1  system, 

P(q)(q  ~  /(«,«))  =  0  (2.29a) 

9(q)  =  0  (2.29b) 

obtained  by  applying  the  coordinate-splitting  matrix  P(q)  to  (2.28a),  is  also  well- 
developed.  By  the  construction  of  P(q),  it  is  easy  to  see  that  the  solution  of  the 
CS  iteration  is  equivalent  to  that  of  the  local  state-space  ODE  of  the  independent 
coordinate  x. 


The  convergence  of  the  CS  iteration  can  be  carried  out  on  a  smooth  constraint 
manifold  M.  Assume  that  for  any  go  €  M,  there  exist  X  €  JRpxn  and  Y  €  2Rmxn 
such  that 

ll(G(9)y)~1||  <  C,  (2.30) 

l|G(9i)T-G(S2fll<C2||S!-92||  (2.31) 

for  some  C\  and  C2,  where  g,  gi,  and  <72  are  in  a  neighborhood  [/(go)  of  go-  Applying 
a  linear  discretization  operator  ph  with  stepsize  h  to  (2.29)  yields  a  nonlinear  system 


F(q)  = 


'  P(q)r(q)  ' 

.  q(q)  J  ’ 


where  the  residual  function  is 


(2.32) 


r(q,t)  =  ph(q)  -  f(q,t).  (2.33) 

Let  { qj }  be  the  solution  obtained  by  applying  the  simplified  Newton  method  to  (2.32), 
and  {qj}  be  the  solutions  obtained  by  applying  Algorithm  2.1.  Convergence  of  the 
CS  iteration  can  be  shown  under  the  conditions  (2.30)  and  (2.31).  Suppose  that 
the  iterative  solutions  {g,}  from  applying  the  simplified  Newton  method  to  (2.32) 
converge  to  g“.  For  go  =  go,  {qj}  generated  by  the  CS  iteration  for  (2.32)  satisfies 

Axk  =  (I  +  Q(q0)TQ(qo))Axk, 


where  Axk  is  defined  by  the  simplified  Newton  iteration.  Consequently,  {g;  }  converges 
to  q“  as  {g;  }  did. 
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Convergence  of  the  CS  iteration  can  be  assured  for  a  sufficiently  accurate  ini¬ 
tial  guess.  Another  sufficient  condition  for  convergence  requires  that  the  numerical 
integration  satisfies 

..(dr  dCf'sY1  ,,  ^ 

•U+Tn  IISCo<1  <2-34) 

for  some  Co  and  s  =  —(GY)~TYTr  in  a  neighborhood  of  q*.  This  implies  an  upper 
bound  on  the  stepsize  h.  For  linear  multistep  integration,  e.g.,  =  the  stepsize 

must  satisfy  9  n 

H  {h  ~  dq^  ~  GTs^J  II  -  1  (2-35) 

where  f3  is  the  leading  coefficient  of  the  numerical  integration  formula. 

For  highly  oscillatory  dynamic  systems  with  a  wide  frequency  band,  the  error 
tolerance  TOL  of  the  numerical  solution  often  satisfies  TOL  <  max^  ||f^||.  Apply¬ 
ing  a  stiffly  stable  numerical  method  to  (2.29),  such  as  BDF  of  order  <  2,  one  may 
take  a  larger  stepsize  to  follow  the  trajectory  of  the  equilibrium,  i.e.,  /  -  GTs  =  0. 
However,  convergence  of  the  Newton  iteration  requires  (2.35),  a  further  restriction 
on  the  stepsize.  Depending  on  how  close  the  predictor  is  to  the  equilibrium  of  highly 
oscillatory  components,  the  Newton  direction  imposed  by  the  Jacobian  can  excite 
the  high-frequency  oscillations.  When  applying  the  Newton  method  directly  to  the 
discretization  of  (2.28),  an  even  more  severe  problem  in  Newton  convergence  is  ob¬ 
served,  and  illustrated  by  the  numerical  experiments  in  Section  5.  The  limitation  on 
the  stepsize  due  to  the  Newton  convergence  failures  for  highly  oscillatory  nonlinear 
multibody  systems  can  be  overcome  via  a  modification  of  the  CS  iteration  which  we 
call  the  CM  iteration. 


3  Highly  Oscillatory  Systems  and  the  CM  Itera¬ 
tion 

3.1  The  CM  iteration 


In  large-scale  multibody  mechanical  systems,  most  of  the  unwanted  oscillations  are 
due  to  the  noise  of  high-frequency  forces,  where  the  amplitude  is  well  below  the 
solution  tolerance.  However,  small  perturbations  in  the  position  can  cause  drastic 
changes  in  the  Newton  direction.  This  results  in  difficulties  for  convergence  of  Newton- 
type  methods.  To  remedy  this  problem  in  the  CS  iteration,  we  reduce  the  noise  from 
the  oscillations  by  setting  =  0  in  the  Newton  iteration  matrix,  since  it  is  the 

main  source  contributing  to  the  rapidly  changing  Newton  direction.  The  basic  idea 
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of  the  CM  iteration  is  to  approximate  the  Newton  direction  of  (2.29)  via  an  oblique 
projection  to  the  unconstrained  ODE 

P(qo)(q~f(q,t))  =  0  (3.1) 

for  a  q0  close  to  the  solution  q,  e.g.,  Giqo^iq)  invertible.  Since  P(q0)  is  no  longer 
varying  with  q,  =  0  is  attained.  When  applying  a  stiffly  stable  numerical  inte¬ 
grator  to  highly  oscillatory  problems,  this  modification,  foi  some  go  close  enough  to 
the  smooth  solution,  overcomes  the  difficulties  in  the  CS  iteration. 

Applying  a  discretization  method  to  (3.1)  coupled  with  constraint  (2.29b)  leads 
to  the  nonlinear  system 

m  =  [  p%(9)  ]  •  <3-2) 

The  corresponding  Lagrange  multiplier  form  of  (3.1)  is 

q  +  GT(q0)\  -  f{q,  t)  -  0  (3.3) 

where 

A  =  (G(q0)Y)-TYT(f(q,t)  -  q). 

A  convergence  result  for  the  modified  CS  iteration,  denoted  by  CM ,  is  given  in  the 
following.  A  detailed  algorithm  for  applying  the  CM  iteration  to  (1.1)  is  presented 
at  the  end  of  this  section. 

We  first  give  an  upper  bound  for  the  difference  between  the  derivative  of  the 
projected  vector  P{q)r{q)  and  the  projected  derivative  P(q )^j^. 


Lemma  3.1  Suppose  conditions  (2.30)  and  (2.31)  hold.  Then 

ll-rlp(«M«)l  -  £  eoCiC2||yrr(,)||  (3.4) 

dq  dq 

in  D{q0 ,  g0)  C  U(qo),  where  D(q0,  g0)  is  the  disc  in  Rn  with  center  q0  and  radius  g0- 


Proof.  The  inequality  is  a  direct  consequence  of  (2.9).  Subtracting  (2.9)  from  P{q)fq 
and  taking  the  norm  of  the  remainder  yields 

\\^irmi)\  -  =  iiPMdGwr(^)"rh||. 

Since  the  row  vectors  of  P(q)  are  p  orthonormal  vectors  in  lRn,  applying  the  Cauchy 
inequality  gives 

in,)*«ll  <  IlfSjpl  <  ,.C2||(Gy)-^r(9)|| 

for  all  q  €  D(qo,  Q0)  C  U(q0).  Condition  (2.30)  implies  the  result  in  (3.4).  □ 


14 


3.2  Convergence  of  the  CM  iteration 

An  estimation  of  the  distance  between  the  solutions  of  (2.32)  and  (3.2)  is  presented 
in  the  following. 


Theorem  3.1  Suppose  conditions  (2.30),  (2.31),  and 

II  (|)  ‘  II  <  Co  <  1  (3.5) 

hold  in  a  neighborhood  of  q *  such  that  {^}  generated  by  the  CS  iteration  converges 
to  q‘.  Choosing  q0  =  q0,  the  sequence  {g*,}  generated  by  the  CM  iteration 

Qk+i  =qk-  J(qo)-1  Fo(qk)  (3.6) 

where  J  =  converges  to  q~ .  Furthermore,  the  distance  between  q*  and  q*  is 
bounded  above  by 

nr  -  rii  <  cai^r  )iinr  -  ®ii  +  nr  -  rf + nr  -  ?<,ii2)  (3.7) 

for  some  moderate  constant  C . 


Proof.  Since  J  is  nonsingular  and  its  components  are  smooth  functions,  using  (2.21) 
and  (2.22)  we  can  write 


1  =. 

Ip  +  QTQ 

0 

-1 

Ip  QT(GY)~'  1 

dr 

-Q 

Im  j 

.  0  (GY)-' 

dq 

■1-1 


for  the  CM  iteration,  provided  that  ^  is  invertible.  By  conditions  (2.30)  and  (2.31), 
we  have 


IP  +  QTQ  0 
-Q  An 


-1 


and 


iP  mr1 1 

0  (GY)-' 


II  ^  A3 

<  c4 


for  some  constants  C3  and  C4.  Thus,  the  contractive  condition  (3.5)  implies  conver¬ 
gence  of  the  CM  iteration. 


To  show  (3.7),  we  first  observe  from  (3.3)  that 

r(q*)  +  G(q0)r{q)  =  0 

where 

G(qo)  —  GT  (qo)(G(qo)Y)~TYT . 


(3.8) 
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Also,  g“  is  a  solution  of  (3.3),  i.e., 


r(?*)  +  a(«W)  =  0-  (3-9) 

Using  (3.8)  and  (3.9),  and  adding  and  subtracting  £?(goM9*)  yields 

(/  +  0(flb))(r(f )  -  rtf))  =  (G(flb)  ~  W)M«T)- 


Applying  the  mean- value  theorem  to  the  differentiable  functions  r  and  Q ,  we  have 

(/  +  C(« o))^(«)(r  -  «*)  =  -  «■) 

for  some  q,  i  =  1,2.  Premultiplying  the  above  equation  by  XT  yields 
^(9o)^(9i)(9*  -  9*)  =  XT^(q2)r(qm){q0  -  qm). 

Since  5(5")  -  g(g*)  -  0,  we  obtain 


\p(qo 

.  G  (qo 


(r  -  o  = 


(«W<r)(»  -  «■) 

§(«o)(9‘  -  5o)2  +  ^ <9o)(f  -  9o)2 


using  the  expansion  of  5(9’“)  and  g(g*)  around  go-  From  (2.21),  (2.22)  and  the  as¬ 
sumption  of  an  invertible  >  we  can  write 


II q  -  9*11  <  C5||AT^r(g')||||g0  -  9*11  +  C6(||g*  -  9o||2  +  II q  ~  9o||2) 
for  some  C5  and  C$.  This  implies  (3.7).  □ 


Note  that  (3.5)  for  the  CM  method  is  analogous  to  (2.34)  for  the  CS  method.  In¬ 
stead  of  (2.35)  for  multistep  integration  methods  using  the  CS  iteration,  the  stepsize 
condition  for  the  CM  iteration  is  given  by 


(3.10) 


in  accordance  with  (3.5). 

It  is  easy  to  see  that  {g*.}  of  the  CS  iteration  and  {g*}  of  the  CM  iteration  are 
the  same  if  the  constraints  g(g)  are  linear.  In  general,  the  rate  of  convergence  of  the 
CM  iteration  is  superlinear ,  using  the  Dennis-More  Characterization  Theorem  [5]. 
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3.3  Implementation  of  the  CM  iteration 


The  CM  iteration  can  be  implemented  via  the  same  procedures  as  those  in  Algorithm 
2.1,  but  the  computational  cost  is  considerably  lower.  First,  the  Jacobian  of  (2.13) 
becomes 


{Qni  t>n) 


dph  Qn 

d(M{qn  )phVn)  _  dfn 

dQn  dqn 


-I 

MdehUa.  _  §Us. 

dvn  dvn 


(3.11) 


since  P(q )  is  held  fixed  in  the  nonlinear  equations.  Comparing  (3.11)  to  the  Jacobian 
(2.16),  notice  that  the  terms  due  to  the  derivatives  of  P(q)r  have  vanished.  Second, 
the  solutions  of  P„A^  and  Pn Auf  at  each  iteration  in  Algorithm  2.1  are  replaced  by 
P(q)4n  and  P(<l)Vn,  which  reduces  the  cost  since  P(q)  has  been  previously  computed. 


Using  the  simplified  Newton  iteration  for  (2.12),  the  CM  iteration  is  carried  out 
as  follows: 


Algorithm  3.1  [CM  iteration] 

Step  1.  Apply  (2.2)  to  GT(q )  and  Ji{q,v)  as  in  (3.11). 

Step  2.  Solve  for  Aq^  and  Av^  in  (2.26). 

Step  3.  Compute  P(q)Aq$  and  P(q)Av„. 

Step  4.-6.  same  as  Algorithm  2.1. 

Remark  3.1  When  applying  the  CM  method  to  (1.1),  the  residual  function 

r2(t\  q ,  t)  =  M(q)phv  -  f(v,  q,  t ) 

should  be  replaced  by 

r2(v,q,t)  =  phv  -  M(q)~1f(v,  q,t) 

for  a  nonsingular  M{q). 


4  Rate  of  convergence  for  highly  oscillatory  multi¬ 
body  systems 


High  frequency  oscillatory  forces  often  appear  in  the  modeling  of  vehicle  suspension 
systems,  modal  analysis  in  structural  dynamics,  or  modeling  oscillations  in  computer- 
aided  engineering  etc.  For  simplicity,  we  consider  the  constrained  dynamic  system  of 
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0 

0 


(4.1a) 

(4.1b) 


(1.1)  with  a  dominant  oscillatory  force 

M(q)v  +  Gt\  +  ^  77(g)  -  f{v,  g,  t )  = 

9(<l)  = 

where  I  may  be,  for  example,  the  coefficients  of  stiff  springs;  i.e.,  0  <  f  <  1.  In 
practice,  77(g)  is  usually  oblique  towards  KerP(q),  i.e.,  the  oscillatory  force(s)  acts  on 
both  the  independent  and  the  dependent  coordinates.  For  the  purpose  of  obtaining  a 
smooth  solution  with  large  stepsizes  solving  those  types  of  problems  [15],  we  will  show 
that  the  CM  iteration  can  be  very  effective  for  many  classes  of  nonlinear  oscillatory 
forces. 

In  the  modeling  of  deformable  multibody  systems,  the  nonlinear  oscillatory  forces 
in  (4.1)  are  usually  derived  from  the  theory  of  linear  elasticity,  i.e.,  for  some  functions 
q  such  that  the  oscillatory  forces  may  be  written  as  Ig.  We  can  use  these  functions 
g  to  write  the  nonlinear  force,  e.g. 

\v(q)  = 


and  then  append 


q  -  77(g)  =  0 


to  the  constraint  equations.  The  oscillatory  forces  then  become  linear  with  respect 
to  the  variables  g.  In  fact,  if  the  oscillatory  forces  are  produced  by  a  finite  element 
approximation  of  the  deformation  of  bodies,  components  of  g  are  associated  with 
some  body-fixed  local  coordinates  via  the  orientation  transformation  matrix,  whose 
entries  often  are  slowly  varying  in  time. 


Deformation  forces  are  the  most  common  potential  forces  that  can  produce  small 
amplitude  high-frequency  oscillations,  and  they  are  usually  linear  writh  respect  to  the 
local  coordinates  [4,  23].  For  these  reasons,  we  consider  the  class  of  oscillatory  forces 
in  the  form 

77(g)  =  B(t)(q  -  b0(t ))  (4.2) 

where  the  components  of  B  and  60  are  slowly  varying.  In  particular,  B  and  bo  may 
be  functions  of  some  constraint-driven  generalized  coordinates.  For  example,  £?(#)  in 
the  2D  bushing  problem  in  [22]  has  the  form 


COS# 

sin# 

0  ‘ 

kx 

0 

0  1 

cos# 

sin# 

0  ’ 

—  sin# 

COS# 

0 

0 

ky 

0 

—  sin# 

cos# 

0 

0 

0 

1 

— ky  sin#  kx  cos  6 

ke 

0 

0 

1 

(4-3) 

where  kx,  ky  and  k°  are  positive  constants.  When  9(t)  is  smooth  or  constrained, 
assumption  (4.2)  is  valid. 
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system'can  eqUati°nS  °f  m°t!°n  °f  the 

M(q)v  +  -B(q  -  b0 )  +  Gt A  -  f(v,  q,t)  =  0  (4.4) 

assume*  that' "  Fy°”  aSSUmption  <2'31)  «"*traint  manifold,  we  can  also 


for  all  q. 


1 

-  >  max 

f  ll«*l  ll,l|l<3  Il=i 


I  dG(q)ui 
1  dq 


«2ll 


(4.5) 


In  the  context  of  the  CS  iteration,  the  problem  of  convergence  of  the  Newton 

by  “g  the  “  ~  £  %% 

V(j)  =  g(q)T(GY)-TYTr  (4  6) 

/  Mij  -B(q  bo).  The  reduced  potential  force  generated  by  (4.6)  is 


=  ^  =  CT(GK)-^r. 


(4.7) 

At  each  iteration,  the  reduced  potential  force  acts  along  the  normal  direction  of  the 
constraint  manifold.  The  gradient  of  the  correction  term  yields 


-r  /  ^/i-rvTi  d^(9)s 


V2V'(q)  =  (/  -  GT(GY)~TYT) 
where  s  =  (GY)~TYTr.  Applying  KT  to  (4.8),  gives 


dq 


(4.8) 


and  applying  XT  to  (4.8)  yields 


YTV2Vr(q)  =  yT(/  _  GT{GY)~TYT)^Ml  _  o 


XTV2Vr(q )  =  P(q)^LM 


dq 

When  high-frequency  oscillations  appear  in  the  system  e  e  e  _»  n  th*  r  ,1  j 
potential  force  also  becomes  oscillatory  if  Y^r  is  nonzero.  This  is  the  general  c^ 
hen  the  solution  is  not  at  an  equilibrium  position.  Nevertheless,  convergence  of  the 
CS  iteration  can  be  achieved  by  using  a  small  enough  stepsize,  e.g.,  hfZy/i. 

Theorem  4.1  Let  (q,v)  be  the  solution  of  the  nonlinear  system  which  results  from 

rar 2kj  uppose  the  ;:^z 

rq  -itLn-r  W  Ml  ’  a  "  °"  ~  and  J(Qo)  is  nonsingular.  Then  the 

CS  iteration  converges  if  h2  <  ce  for  some  moderate  c 
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Proof.  For  the  convergence  of  the  CS  iteration,  we  need  to  show  that  (2.34)  >s  valid, 
where  r(j)  is  defined  in  (4.4).  For  (2.34),  we  have 


=  ll^rr1  -  hi +  °(h) 


where  /?  >  0  is  the  leading  coefficient  of  ps,  and  ||M||  is  not  zero.  Consequently,  for 
€  <C  1,  (2.34)  is  valid  provided  that  h  «  y/e.  □ 

Under  the  conditions  of  the  above  theorem,  a  convergence  result  for  the  CM  itera- 
Uon  can  be  obtained  provided  the  assumptions  of  Theorem  3.1  are  valid  In  many 
applications,  following  the  oscillations  is  not  of  interest,  taste  ad, , 
laree  time  step  to  damp  out  the  oscillations  of  small  amplitude  but  high  frequen  y. 
For  this  reason,  we  now  consider  only  the  multistep  numerical  integration  methods 
that  are  strictly  stable  at  infinity  and  4-stable,  such  as  the  lower- order  (.  .e  2 

BDF  methods  [11]  The  convergence  of  L-stable  implicit  Runge-Kutta  methods  to  the 
“™f  the  highly  oscillatory  ODE  of  multibody  mechanical  systems  can 
bTfoundt  M  Here  we  focus  on  the  convergence  of  the  CM  iteration  for  constr  ained 
multibody  systems  with  oscillatory  forces  when  applying  the  above-mentioned  line 
multistep  methods. 

Numerical  solutions  on  the  slow  manifold  can  be  evaluated  using  the  equilibrium 
of  (4.1),  i.e. ,  the  slow  solution  [2,  13]  satisfies 

n(q)  -  €(/(«,  q)  -  Wr(g)  -  M(q)v)  =  0, 

and  the  smooth  solution  is  its  asymptotic  expansion  to  some  order  of  e  mound  the 
manifold  { q  I  n(o)  =  0).  In  the  linear  form,  the  smooth  solution  of  (4.1)  not 
™”om  B(,  -  6„)  =  o’ since  }  >  &  For  the  strongly  damped  numerical 
R/  _  h  }  Ole)  as  t  — i  oo  During  the  iterative  solution  onto  the 
itTlnifomX  constraints  may”  not  be  satisfied,  which^s  ^-action 

b* the  red7td0 

lotentfal  This  yields  a  superior  performance  of  the  CM 

the  CS  iteration  for  computing  the  smooth  solution  of  (4.1).  The  result  is  explained 
in  the  following. 

Lemma  4.1  Let(f,V)  be  the  smooth  solution  of  (U)*')  linear  ana Ihthe ■step, 
size  of  the  multistep  integration  method.  Suppose  the  startmg  values  (go,  Vo)  for  («  ,v 
on the  smooth  solution  of  (U),  i.e.,  M  =  0(e)  and  r  (?,*)  =  0(h),  sattsfy  (2.30), 

(2.31)  and  M(so)ps(uo)  -  /(fo,9o)  =  °(h)  <4  9’' 
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where  ph  is  the  corresponding  discretization  operator.  Applying  the  CS  and  CM 
iterations  to  (4-1)’  approximate  Jacobian  matrix  for  the  CS  iteration  satisfies 

l|J(«>,»»)  -  ■/(«-,  f')ll  =  6-0(h)  +  0(h)  (4.10) 

where  <5  =  \\Bqo  —  Bq” ||,  and  J(q,v )  is  the  Jacobian  of  (4-1)-  For  the  CM  iteration, 
we  have 

lk(?o,^o)  ~  ||  =  0(h)  (4.11) 

where  J  is  the  approximate  Jacobian  in  the  CM  iteration. 

Proof.  The  difference  between  the  Jacobian  at  (qo,vo)  and  (q“,vm)  can  be  written  as 

\\Jo  -  J-||  <  ||P(?o)^(?o)  -  P(g-)^q( ?-)ll  +  II^(?oM9o)  -  ^(s'W?-)ll  +  0(h) 

since  the  initial  values  satisfy  (4.9).  Under  the  conditions  (2.30)  and  (2.31),  we  may 
choose  common  X  and  Y  for  P(qo)  and  P(q’)  such  that  the  first  term  on  the  right- 
hand  side  of  the  above  inequality  can  be  rewritten  as 

l|P(9)(|(?o)  -  ^(?-))ll  =  0(h) 

for  some  q  €  [<?0)9“]i  since  =  \B  +  0(h)  allowing  the  cancellation  of  -(B.  The 
second  term  yields 

ll^(9o)r(<Jo)  -  ^(9-)r(9‘)||  <  ||r(g„)  -  r(,-)||0(h)  =  l||B9o  -  Bq'\\0(h) 

according  to  Lemma  3.1.  Thus,  (4.10)  is  proved.  Recalling  J(qo,Vo)  from  (3.11),  we 
have 

l|.Wl<l|f>llll0ll  =  O(ft). 

using  again  Lemma  3.1.  □ 


Theorem  4.2  For  the  initial  values  (qo,vo),  suppose  the  conditions  in  Lemma  4-1 
hold.  Suppose  that  the  CS  and  the  CM  iterations  are  carried  out  by  applying  a 
simplified  Newton  method,  where  the  iteration  matrix  is  computed  at  the  starting 
values  (qo,vo).  If  both  iterations  converge,  then  the  rate  of  convergence  of  the  CS 
iteration  o(CS^  compared  to  that  of  the  CM  iteration  a(CM)  is  given  by 

=  -0(h)  +  0(h) 

where  6  =  || B(q0  —  q*)\\,  and 

=  0(h). 
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Proof.  Since  we  apply  the  simplified  Newton  iteration,  the  solution  of  the  CS  itera¬ 
tion  can  be  written  as 


Qk+ 1 
Vk+l 


H(qk,v  k) 


where 


H(q,  v)  = 


Similarly,  the  CM  iteration  can  be  written  as  the  fixed-point  iteration  of  the  function 


H(q,v)=  l  -J^F(q,v). 


Applying  the  Contractive  Mapping  Theorem,  see  [5]  pp.  93-94,  we  obtain  the  rates 
of  convergence  of  the  CS  and  CM  iterations 


Ccs  =  11/  -  t>-)||  =  IIJo-’M,  -  Oil 


and 

ccm  =  ii/  -  on  =  n/0-'(/o  -  on 

where  J(qm,v’')  =  J“  is  the  Jacobian  at  the  solution  of  the  discretized  system,  and 
the  superscripts  denote  the  respective  iterations.  For  J^1,  we  have 


Jo 


l  +  0(h)  -i 
f B  +  0(h)  |jV/  +  0(h) 


When  e  — ►  0,  the  dominant  components  of  J0  1  are  of  0(1).  From  Lemma  4.1,  the 
rates  are 

aC5  =  - 0{h )  +  0(h) 

and 

acs  =  0(h), 

since  Jq 1  has  no  component  of  0(|).  □ 


5  Numerical  Experiments 

5.1  Point-mass  with  oscillatory  force 

The  first  example  is  a  simple  constrained  multibody  system  under  the  influence  of 
a  highly  oscillatory  force.  Consider  a  unit  point-mass  constrained  to  the  2D  unit 
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circle,  using  q  =  [x,y]T,  the  velocity  v  =  ^  =  [w,z]T,  and  the  constraint  equation 
g(q)  =  +  y2  —  1).  The  equations  of  motion  are 

w  +  x\  —  fx  =  0 
z  +  yX-  fy  =  0 

where  /  =  [/x,/y]T  is  the  applied  force.  Differentiating  the  constraint  g{q)  =  |(x2  + 
y 2  —  1)  twice  with  respect  to  time,  an  explicit  form  of  the  multiplier  A  is  obtained  by 

x  =  -^-2(xfx  +  yfy  +  w2+z2)- 

For  a  highly  oscillatory  force  /(g),  one  can  see  that  X(t)  is  oscillating  with  the  fre¬ 
quency  of  f(q),  and  with  amplitude  proportional  to  the  magnitude  of  ||/||. 

The  numerical  experiments  are  carried  out  using  BDF  of  order  <  2  in  DASSL  [3], 
where  the  local  error  estimation  has  been  modified.  For  the  stabilized  index-2  DAE 
(1.3)  denoted  by  GGL ,  the  local  error  is  estimated  using  only  the  position,  i.e.,  q. 
Moreover,  we  have  also  included  some  experiments  where  the  Newton  convergence  test 
of  GGL  has  been  modified  to  exclude  the  multipliers.  The  corresponding  numerical 
solution  is  denoted  by  GGL*.  For  the  coordinate-split  and  modified  coordinate-split 
iterations,  denoted  by  CS  and  CM,  respectively,  the  local  error  is  estimated  using 
the  independent  variable  XTq ,  as  recommended  in  Section  2.3.  The  CM  iteration 
updates  the  matrix  P(q)  when  a  new  Jacobian  is  required. 

Linear  oscillation 


Let  a  unit  gravitational  force  act  along  the  negative  y-direction,  and  apply  a  linear 
oscillatory  force 

/  = 


l-x 

€ 


L  € 


My+  1)  -  1.0 


There  is  a  stable  equilibrium  at  q  =  [0,  -1]  and  v  =  [0, 0].  The  natural  frequency  of 
the  system  is  u>  =  and  no  dissipative  force  is  present. 


The  numerical  solution  has  been  carried  out  with  a  moderate  solution  tolerance 
ATOL  =  RTOL  =  TOL  =  10~3.  For  a  0  to  0.25  second  simulation,  the  results  of 
several  combinations  of  the  stiffness  coefficient  e  are  presented  in  Table  5.1,  where 
the  initial  values  are  q  =  [0.0, —1.0]  and  v  =  [1.0,0].  The  CM  iterations  show  better 
efficiency  than  those  of  CS,  GGL  and  GGL *  in  all  cases,  i.e.,  comparing  the  numbers 
of  function  and  Jacobian  evaluations  in  Table  5.1.  In  the  Table,  etfs  and  ctfs  denote 
the  number  of  failures  of  the  error  test  and  Newton  convergence  test,  respectively,  in 
DASSL.  Comparing  the  results  of  GGL*  with  those  of  GGL ,  we  observe  an  improved 
Newton  convergence.  As  e  -»  0,  i.e.,  for  higher  frequencies  of  the  oscillation,  the 
CM  iteration  becomes  even  more  efficient.  In  Figure  5.1,  we  plot  the  total  energy  of 
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Table  5.1:  Results  of  the  Constrained  Point-Mass  with  a  Linear  Oscillatory  Force 


Figure  5.1:  Total  Energy  Comparison  of  Linear  Oscillatory  Force  Example;  e  =  10  6 
and  TOL  =  10-3 
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Table  5.2:  Results  of  the  Constrained  Point-Mass  with  an  Oscillatory  Linear  Spring 
Force 


each  numerical  solution.  The  CM  iteration  achieves  the  strongest  damping  because 
DASSL  is  able  to  increase  the  stepsize  faster  with  the  CM  iteration. 

Linear  spring  force 

In  the  next  test,  we  replace  the  linear  oscillatory  force  in  the  previous  constrained 
system  by  a  spring  force 

1  l  ~  x  —  Xo 

e2  l  V  ~  2/o  . 

where  l  =  for  (xo,  Vo)  the  attachment  point  of  the  spring,  Iq  the 

natural  length,  and  j?  the  stiffness  coefficient,  as  shown  schematically  in  Figure  5.2. 
For  unit  mass  and  unit  gravitational  force,  we  set  the  spring  attached  to  (:ro,  j/o)  — 
(0,  —0.5)  and  the  natural  length  lo  =  0.4,  such  that  the  equilibrium  is  at  (0,  —1,0, 0). 

Using  the  initial  conditions  [0.04471,  -0.999, 0, 0],  the  results  of  the  0-0.05  second 
simulation  by  the  GGL,  GGL*,  CS,  and  CM  iterations  are  shown  in  Table  5.2. 
Because  DASSL  is  able  to  increase  the  stepsize,  and  hence  damp  the  solution  faster 
with  the  CM  iteration,  the  CM  method  is  quite  effective  in  these  tests.  In  Figure 
5.3,  we  plot  the  total  energy  of  each  solution.  The  numerical  solutions  of  x,  w,  and  X 
are  presented  in  Figure  5.4. 
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Figure  5.2:  Constrained  Point-Mass  with  a  Linear  Spring 


Total  Energy  v.r  Time 


Figure  5.3:  Total  Energy  Comparison  of  Oscillatory  Spring  Example;  e  =  10  4  and 
TOL  =  10-4 
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5.2  Two-body  pendulum  with  bushing 

The  second  example  is  a  two-body  pendulum  in  2D  Cartesian  coordinates.  Six  gen¬ 
eralized  coordinates,  q  =  [x\, yi,9i,X2, 92,62V  •>  locate  the  centers  of  mass  and  the 
orientation  of  the  bodies.  The  first  body  is  grounded,  and  the  second  body  is  con¬ 
strained  such  that  the  distance  between  a  point  A  of  the  first  body  and  another  point 
B  of  the  second  body  is  fixed,  and  its  orientation  is  held  constant.  This  leads  to  five 
constraint  equations 


9\  =  Xy 

(5.1a) 

92  =  2/1 

(5.1b) 

93  —  $1 

(5.1c) 

g4  =  dABTdAB-lAB 

(5.1d) 

1 

<N 

II 

(5.1c) 

where  lAB  and  9  are  constant,  and 

jab  _  [  £ i  +  Qi  cos  6>i  —  02  sin  6\  —  (x^  4-  &i  cos  62  —  62  sin  #2) 

—  y\  +  Gi  sin#i  +  02 cos 9\  —  (2/2  +  61  sin 62  +  &2COS02)  . 

such  that  A  =  [oi,o2]  and  B  =  [61,62]  in  the  local  reference  coordinate  systems  of 
body  1  and  body  2,  respectively.  In  this  example,  we  use  A  =  [0,0],  B  =  [0,0], 
lAB  =  1,  and  9  =  0  for  the  constraint  equations  (5.1). 

We  apply  a  nonlinear  oscillatory  force  formulated  using  nonlinear  beam  theory  [4]. 
This  type  of  force  arises  commonly  in  flexible  multibody  dynamics  [23].  As  described 
in  (4.2),  the  deformation  force  between  the  i th  and  j th  components  is  a  function  of 
the  relative  displacement  of  the  reference  frames  X\-Y[-Z\  and  Xj-Yj-Zj,  as  shown 
schematically  in  Figure  5.5.  Typically,  the  relative  displacement  is  measured  by 

d{j  =  r  j  +  Ajs'j  —  r,-  —  A,st  (5-2) 

where  s'  and  s)  are  constant  vectors  to  the  origins  of  the  force  reference  frames  in 
their  respective  local  coordinate  systems,  i.e.,  Xi-Yi-Zi  and  Xj-Yj-Zj ,  where  rt,  r,- 
are  the  corresponding  origins  in  a  global  coordinate  system  and  A,  and  Aj  are  the 
transformation  matrices  from  the  global  to  the  local  coordinate  system  [8].  The 
relative  angles,  ©ti  =  [ipxj,9lj,<t>lJ}T ,  are  calculated  as 

Aij  =  {AiBi)T  AjBj 

ipij  =  ~Aij(  2,3) 

9ij  =  Ay  (1,3) 

*>=arctan(l^2jj 
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where  Aij(k,l )  is  the  component  of  the  k th  row  and  Zth  column  of  A{j.  The  matrix 
Aij  is  the  relative  orientation  matrix  of  two  force  reference  frames,  i.e.,  Bt  and  Bj 
are  constant.  The  relative  velocity  is  the  time  derivative  of  the  relative  displacement 
dij  =  frdij.  and  the  relative  angular  velocity  is  Uij  =  uij  —  where  u and  uij  are  the 
angular  velocities  of  bodies  i  and  j  respectively. 


Figure  5.5:  Deformation  Force  of  a  Flexible  Body 

Using  the  above  defined  notation,  the  force  acting  between  the  ith  and  j  th  com¬ 
ponents  due  to  the  deformation  can  be  written  as 

fa  =  AiBi(K*(AiBi)T  dij  +  C^AiBifdij) 

where  K*  is  a  3  x  3  structural  stiffness  matrix  and  is  the  3x3  damping  coefficient 
matrix.  Similarly,  the  torque  acting  between  the  components  is 

Tij  =  AiBi(KTQij  +  C^AiBifuij) 

where  KT  and  CT  are  analogous  to  Kf  and  Cf.  Note  that  the  force  and  torque  in 
this  form  are  linear  functions  of  the  relative  displacement  (dtJ  ,  0tJ)  and  the  relative 
velocity  (d0-, 


Here,  the  2D  bushing  force  has  stiffness  matrix 


0  0  1 
ky  0 
0  ke 
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where  kx ,  ky,  and  ke  are  0(1),  and  damping  matrix 

[c1  0  0  1 

cf  =  0  cy  0 

0  0  c6  _ 

where  .cf ,  cv,  and  ce  are  0(1).  For  this  2D  bushing  example,  (4.3)  becomes 

1  [  1  0  0  1  [  cos 2  6  cos  6  sin  0  0  " 

B(6)  =  -  0  1  0  +  \kx-ky\  cos  6  sin  6  sin 2  6  0  • 

■  €  L  0  0  k6  \  [  0  0  0. 

The  attachment  points  of  the  force  device  are  Sj  =  [0.5, 0]  and  s'2  —  [—0.5, 0]  in  the 
body-fixed  reference  frames  of  bodies  1  and  2,  respectively.  The  bushing  force  intro¬ 
duces  oscillatory  applied  forces,  causing  small  oscillations  of  the  numerical  solution, 
and  yielding  highly  oscillatory  multipliers  in  the  index-2  DAE  (1.3).  The  multipliers 
associated  with  the  highly  oscillatory  components  exhibit  high-frequency  oscillations 
with  large  amplitude.  The  standard  convergence  test  of  the  Newton  iteration  depends 
heavily  on  these  multipliers.  Therefore,  we  modified  the  convergence  test  in  DASSL 
to  exclude  the  test  for  the  multipliers.  In  addition,  the  multipliers  are  computed 
in  the  GGLT  by  applying  the  pseudo-inverse  (GY)~TYT  to  r\  =  —  PhQ^  and  to 

r2  =  v^°\t)—MphV^°\  where  (^°\  t/°))  is  the  predictor  in  DASSL,  and  ph  is  the 

discretization  operator  of  BDF.  The  local  error  is  estimated  by  the  predictor-corrector 
difference  of  {XT q,  XTv)  for  C5,  CM,  and  GGL ',  and  of  ( q ,  v)  for  GGL. 

Using  the  initial  values  q  =  [0,0,0,9.9989e-l,-1.4852e-2,0]  and  v  =  [0,0,0,-6.75e- 
5,-4.5444e-3,0],  numerical  results  with  ATOL  =  RTOL  =  TOL  are  shown  in  Table 
5.3,  where  e  =  10-3,  kx  =  ky  =  ke  =  1,  and  cx  —  cy  =  ce  —  10.  For  this  moderate 
stiffness  e  =  10-3,  all  the  methods  perform  well.  Note  that  the  constraint  violation 
with  these  initial  values  is  O(10-3).  This  implies  that  the  difference  of  the  constraint 
reaction  force  at  the  initial  value  qo  and  that  of  a  q*  on  the  constraint  manifold  is 
<5  «  ||5o  -  9*||  =  O(10-2),  e.g.,  6  =  O(v/10-3).  According  to  Theorem  4.2,  the  rate  of 
convergence  of  the  CS  iteration  is  proportional  to  \0(h).  Increasing  numbers  of  the 
convergence  test  failures  in  DASSL  are  expected  as  e  — ►  0.  In  this  example,  frequent 
convergence  test  failures  occured  when  e  <  10-5  and  TOL  <  10-3.  We  observe  the 
same  difficulties  in  the  Newton  convergence  of  the  GGL  and  GGL *  iterations.  On 
the  other  hand,  the  CM  iteration  with  its  better  Newton  convergence,  as  explained 
by  Theorem  4.2,  is  able  to  take  much  larger  time  steps  and  the  nonlinear  oscillation 
is  damped  effectively.  In  Table  5.4,  the  results  of  e  =  10-6  ,  kx  =  ky  =  ke  =  1,  and 
(f  =  cy  =  c°  =  10  are  shown.  In  Figure  5.6,  we  plot  the  stepsize  taken  by  DASSL  for 
GGL ,  GGL *,  CS ,  and  CM  using  the  stiffness  coefficient  c  =  10-5. 
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